# scottknott.auk
# Tim Menzies
# Sept 12, 2012

"""

stuff

"""

@include "lib.awk"
@include "readnums.awk"
@include "bootstrap.awk"

# small 0.17 medium 0.6 large 1.4 }

function _scottknott(   o) {
  args("-d,data/gots/sc1,-f,%3s,-w,30,-l,0,-h,100,-s,0.17",o) #0.38
  sk(o["-d"],o["-f"],o["-w"],
     o["-l"],o["-h"],o["-s"])
}
function sk(f,   fmt,width,lo,hi,small,    datas,tmp) {
  readnums("cat " f "/*",datas) 
  scottknott(datas,small)
  tiles(datas, fmt,width,lo,hi)
}
function tiles(datas, fmt,width,lo,hi,\
               data,chops) {
  pairs("0.1,-,0.3, ,0.5, ,0.7,-,0.9, "          ,chops)
  tiles1(datas,0,chops,fmt,width,lo,hi)
  for(data in datas) 
    if (data != 0)
      tiles1(datas,data,chops,fmt,width,lo,hi)
}
function overallSd(data,  a) {
  allData(data,a)
  return l2sd(a)
}
function allData(data,a,   n,i,j) {
 for(i in data)
    for(j in data[i]["="])
      a[++n] = data[i]["="][j]
}
function tiles1(datas,data,chops,fmt,width,lo,hi,   n,i) {
  print "\n----| " data " |----------------------"
  n = length(datas[data])
  for(i=1;i<=n;i++)  {
    print sprintf("%3s,%20s,",
                  datas[data][i]["rank"],
                  datas[data][i]["rx"]),
      tile(datas[data][i]["="],
           chops,width,fmt,lo,hi)
}}
function scottknott(datas,small,  data) {
  for(data in datas) 
    sk1(datas[data],small,1,1,length(datas[data]),0) 
}
function sk1(rxs,small,rank,lo,hi,  cut,i) { 
  cut = sk0(rxs,small,lo,hi)
  if (cut) {
    rank = sk1(rxs,small,rank,lo,cut-1) + 1
    rank = sk1(rxs,small,rank,cut,hi) 
  } else {
    for(i=lo;i<=hi;i++)
      rxs[i]["rank"] = rank
  }
  return rank
}
function sk0(a,small,lo,hi,\
             s,  s2, n, mu,                     \
             ls,ls2,ln,lmu,                     \
             rs,rs2,rn,rmu,                     \
             i,best,tmp,cut) {
  for(i=lo;i<=hi;i++) {
    s += a[i]["s"]; s2+= a[i]["s2"]; n += a[i]["n"]
  }
  mu= s/n
  best = -1; 
  for(i=lo+1;i<=hi;i++) {
    ls  += a[i-1]["s"]; 
    ls2 += a[i-1]["s2"]; 
    ln  += a[i-1]["n"]
    rs   = s - ls; 
    rs2  = s2- ls2; 
    rn   = n - ln; 
    rmu  = rs/rn
    lmu  = ls/ln
    tmp = ln/n * (mu - lmu)^2 + rn/n * (mu - rmu)^2 
    if (tmp > best) 
      if  (ghedge(rmu,lmu,rn,ln,rs,ls,rs2,ls2) > small) 
        if (significantlyDifferent(a,lo,i,hi)) {
            best = tmp; cut = i }}
  return cut
}
function ghedge(rmu,lmu,rn,ln,rs,ls,rs2,ls2,\
              n, ln1,rn1,rsd, lsd,sp,correct,g) {
  n       = rn + ln
  ln1     = ln - 1
  rn1     = rn - 1  
  rsd     = sd(rn,rs,rs2)
  lsd     = sd(ln,ls,ls2)
  sp      = sqrt((ln1*lsd^2 + rn1*rsd^2)/(ln1 + rn1))
  correct = 1 - 3/(4*(n - 2) - 1)
  return abs(lmu - rmu)/sp * correct
}
function significantlyDifferent(a,lo,cut,hi,\
                               i,j,x,r,l,rhs,lhs) {
  for(i=lo;i<=hi;i++) 
    for(j in  a[i]["="]) {
      x = a[i]["="][j]
      i < cut ?  lhs[++l] = x :  rhs[++r] = x
    }
  print ":lo",lo,":cut",cut,":hi",hi,":left",l,":right",r > "/dev/stderr"
  return bootstrap(lhs,rhs,500,0.05)
}
